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Abstract 

This paper presents the results of a numerical study of the direct contact condensation on a slowly moving hori- 
zontal liquid surface. The geometrical configuration and the input conditions used to obtain numerical solutions 
are representative to those of experiments of Celata et al. 7 The effects of Prandtl number (Pr), inflow Reynolds 
number (Re in ), and Richardson number (Ri) on the condensation rate are investigated. Numerical predictions of 
condensation rate for laminar flow are in good agreement with experimental data. The effect of buoyancy on the 
condensation rate is characterized by Richardson number. A correlation based on the numerical solutions is devel- 
oped to predict the average condensation Nusselt number in terms of Richardson number, Peclet number, and 
inflow Reynolds number. 


Nomenclature 

A surface area 

Ar aspect ratio, H/D 

B tank to jet diameter ratio, D/d 

C p specific heat at constant pressure 

d,D jet diameter, tank diameter 

g gravitational acceleration 

H liquid height 

h c heat transfer coefficient at the interface h c * k( 3T/3xy(T s - T in ) 

h^ latent heat of condensation 

Ja Jakob number, C p (T s - T^J/h^ 

k thermal conductivity 

m c local condensation mass flux 

Nu c Condensation Nusselt number, h c D/k 

Nu co local condensation Nusselt number for Ri * 0 
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p,p g pressure, equilibrium hydrostatic pressure 

p* dimensionless pressure, (p - Pg)/Puj n 2 

Pe Peclet number, Pe = Re^Pr 

P Prandtl number, Pr = C p p/k 
inflow volume flow rate 

r radial coordinate measured from the centerline 
r* dimensionless radial coordinate, r/D 

Re^ inflow Reynolds number, pu^ (D - 2d )/p, 

Ri Richardson number, Ri * P (T s - T in )gH/u eq 2 
T temperature 

T* dimensionless temperature, (T - T« n )/(T S - T^) 
u axial velocity 

u c condensation-induced velocity 

u^ equivalent horizontal velocity defined in Eq. (9) 
u* dimensionless axial velocity, u/u in 

v radial velocity 

v* dimensionless radial velocity, v/u in 

x axial coordinate measured from the tank bottom 

x* dimensionless axial coordinate, x/D 

Greek symbols: 

P thermal expansion 

p. dynamic viscosity 

p liquid density 

Subscripts: 

in . evaluated at inflow location 

out evaluated at outflow location 
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s evaluated at liquid-vapor interface 

c evaluated at condensation condition 


Introduction 


The pressure control of cryogenic tanks during the on-orbit storage and transfer of cryogenic liquids, such as 
hydrogen and nitrogen under microgravity conditions is a challenging problem. A preferred method of tank pres- 
sure control in space environment is the use of a thermodynamic vent system along with fluid mixing. 1,2 Direct 
contact condensation of vapor on the liquid surface plays a central role in such a pressure control system. 

The authors 3 " 5 numerically obtained the condensation rate as a function of the associated system and flow 
parameters in a laminar axial jet-induced mixing tank. One important conclusion of these investigations is that the 
interface condensation rate could be simply determined by the jet volume flow rate and jet subcooling 3,4 for con- 
ditions where buoyancy effects are negligible. The effects of buoyancy on the condensation rate can be character- 
ized in terms of Richardson number. 5,6 Numerical predictions 5 show that in a laminar jet induced mixing system 
the condensation rate is significantly reduced for higher values of Richardson number. The authors 3 ' 5 are not 
aware of any experimental data of the vapor condensation rate on a liquid surface in a laminar jet induced mixing 
system. Therefore, the numerical predictions could not be compared with experiments. 

Recently, Celata et al. 7 conducted an experiment and measured the condensation rate of saturated steam on a 
horizontal, subcooled, slowly moving water surface. The range of liquid flow rates covered in the experiment 
encompasses the laminar and the turbulent flow regimes. A theoretical model that contains three empirical con- 
stants was developed to predict the interfacial heat transfer coefficient. However, the description of the transport 
phenomena in the near- interface region was not presented in the paper. Also, the effects of fluid properties, buoy- 
ancy, and mass flow rate that can generally be expressed in terms of Prandtl number, Richardson number and 
Reynolds number were not shown explicitly. The objective of this study is to obtain a numerical solution of 
Celata's experiment, clarify the effects of the pertinent dimensionless parameters on the condensation rate, and 
compare the prediction with the measurements 7 A brief description of Celata's experiment and its adaptation to 
the present numerical simulation is presented in the following sections. 

The Physical Problem 

The physical system and the coordinates used to analyze the problem are similar to those of Ref. 7 and are 
shown in Fig. 1. A circular cylindrical tank of diameter D contains liquid with a filling height H. A uniform 
annular inflow of velocity u in is introduced into the tank from the outer portion of the tank bottom toward the 
liquid-vapor interface. Liquid is simultaneously withdrawn from the central portion of the tank at the same vol- 
ume flow rate Q in such that the liquid fill level is kept constant. The circular liquid outflow has a diameter, d. 
The tank pressure is maintained constant resulting in a constant interface saturation temperature T s . The inflow 
liquid is subcooled and has a constant and uniform temperature T in . The outflow region is assumed to have zero 
temperature gradient. 


The energy balance at the interface yields the following relation between condensation mass flux and heat trans- 
fer at the interface 


m c hfgA s — 


m c h fg dA s 




dA c 


where m c is the average mass condensation flux over the interface area A s . 
transfer coefficient (h c ) and average condensation Nusselt number (Nu c ) over 

r _ m c h fg 

C = (T, - T m ) 


The average condensation heat 
the interface can be defined as 


( 1 ) 


( 2 ) 


and 
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Mathematical Formulation 


The steady-state, incompressible Navier-Stokes equations are employed to solve the present problem. The grav- 
ity is acting in the vertical negative-x direction and the effect of buoyancy is accounted for by using the Bous- 
sinesq approximation. The dimensionless form of the governing equations are expressed as: 


du* 2 du*r*v* 
+ 

dx* r*dr* 


dp * 
dx* 


du* dr*v* 
dx* r*dr* 


(1 - 2/B) 

a 2 u* + a fdu*) 

B 2 1 - (2/B) 2 

-L_ L 

Rein 

dx* 2 r*dr*^ dr*J 

16Ar 3 


(4) 


RiT* (5) 


du*v* dr*v* 2 
dx* r*dr* 


du*T* 

dx* 

where the parameter B, aspect ratio 
number Ri are defined as 


dp* (1-2/B) v* (1-2/B) 


dr* 


Re;, 


Re;. 


d 2 v* 

3 f 

'r^i 

dx* 2 

r*dr*l 

. dr *1 


9r*v*T* _ (1 - 2/B) 

r*dr* 


Re^Pr 


3 2 T* 




dx* z r*3r*v dr* 


Ar, Prandtl number Pr, inflow Reynolds number Re- , and Richardson 


( 6 ) 

(7) 


B = D/d, Ar = H/D, Pr = C p p/k, Re^ = p Uin (D - 2d)/fl, Ri = (3(T S - T in )gH/ Ue<1 2 (8) 

Other variables and dimensionless parameters are defined in the nomenclature. It is to be noted that in the defini- 
tion of Richardson number (Ri), u^ is defined as an equivalent horizontal velocity with which the flow passing 
through the opening between the interface and the outflow location has the same volume flow rate as subcooled 
inflow, i.e., 

re dHu cq = (9) 

The boundary conditions are required to solve Eqs. (4) to (7). At the centerline, the symmetric conditions are 
used: 


v* = 0, 


du* dT* 
dr* dr* 

At the solid walls, the no-slip and adiabatic conditions are employed: 


= 0 at r* = 0 


( 10 ) 


The free surface is assumed to be flat 
face, u c , is equal to 


u* = v* = 0, 


dT* dT* 


= 0 at walls. 


dx* dr* 

(wave free) and shear free. The condensation-induced velocity at the inter 


( 11 ) 


u c = m c /p (12) 

and is generally negligibly small compared with surface velocity. 3 The temperature of the interface is kept con- 
stant at T s , Thus, the boundary conditions applied at the interface are. 


u* = 0, — = 0, T* = 1 at x* = H/D (13) 

dx* 

For the inflow, the velocity and temperature are assumed to be uniform: 

u* - 1, v* = 0, T* = 0 at the inflow location. (14) 

In order to maintain liquid level constant, the liquid withdrawn at the central nozzle should have a volume flow 
rate equal to the sum of inflow volume flow rate Q in plus induced interfacial condensation flow rate. Since u c 
is assumed negligibly small, the outflow volume flow rate is approximately equal to the inflow volume flow rate. 



Assuming the outflow velocity is uniform, the outflow velocity is obtained by Q[JA mV At the outflow location 
a zero temperature gradient is assumed. Thus, 

dT* 

u* = “Ain/Aout> v* = 0, — - = 0 at outflow locadoa (15) 

dx* 

Based on the above dimensionless variables, the average heat transfer coefficient h c becomes 


h c 


r 0.5 


- 8(k/D) 



and the average condensation Nusselt number is 



(16) 


(17) 


Numerical Method of Solutions 

Equations (4) to (7) are numerically solved by a finite-difference method. The finite-difference equations are 
derived by integrating the differential equations over an elementary control volume surrounding a grid node appro- 
priate for each dependent variable. 8 A staggered grid system is used such that the scalar properties, p* and T*, 
are stored midway between the u* and v* velocity grid nodes. The bounded skew hybrid differencing (BSHD) 
is incorporated for the convective terms 8 and the integrated source terms are linearized. Pressures are obtained 
from a predictor-corrector procedure of the Pressure Implicit Split Operator (PISO) method 9 which yields the pres- 
sure change needed to procure velocity changes to satisfy mass continuity. The governing finite-difference equa- 
tions are solved iteratively by the ADI method with under relaxation until the solutions are converged. 

Calculations are performed with a nonuniform grid distribution with concentration of the grid nodes in the 
centerline, near-wall, and near-interface regions where the gradients of flow properties are expected to be large. 
The nonuniform grid distribution is generated by using an exponential function of Roberts’ transformation. 10 
Several grid distributions have been tried and the grid pattern of 48 by 41 has shown to give reasonably grid- 
independent solutions of interface quantities. Calculations are performed on a CRAY-XMP computer located at 
NASA Lewis Research Center. The convergent solutions are considered to be reached when the maximum of 
absolute residual sums for u*, v*, and T* variables is less than 10' 5 . 

Results and Analysis 

Based on the system sketched in Fig. 1 and the associated governing equations and boundary conditions, the 
average condensation Nusselt number can be expressed as 

Nuc = f(B, At, Ajn/A^,, Pr, Re^, and Ri) (18) 

The geometrical configuration and the input conditions used to obtain numerical solutions are similar to those of 
experiments of Celata et al. 7 such that quantitative comparison between the numerical solutions and experimental 
data can be made. Accordingly, the system parameters (B, Ar, and A ij JA m ) and flow parameters (Re in , Pr, and 
Ri) are specified as 

B = 1333, Ar = 0.075, AJA oxxi = 174 (19) 

and 

36 < Re^ < 330, 1.1 < Pr < 7.1, 0 < Ri < 5 (20) 

Figure 2 shows that for Ri = 0, the average condensation Nusselt number Nu c is increasing with inflow 
Reynolds number (Re in ) and Prandtl number (Pr). A correlation based on the numerical solutions (within ±1 per- 
cent) to predict Nu c at Ri ■ 0 is expressed as: 

Nuco = 3.27(Pe/ 2 (21) 

where Pe is the Peclet number and is defined as Pe = Re in Pr. 


5 



The buoyancy effect is characterized by Richardson number, Ri. As Ri increases, the buoyancy impedes the 
subcooled liquid inflow to reach the interface, resulting in a lower condensation rate. The average condensation 
Nusselt number (Nu c ) as a function of Richardson number for various Pec let numbers and inflow Reynolds num- 
bers are shown in Figs. 3 and 4, respectively, in which nH C 0 is the average condensation Nusselt number with 
Ri equal to zero (no buoyancy effect). For a given inflow Reynolds number and Peclet number, the condensation 
rate is linearly decreasing with increasing Richardson number. The effect of buoyancy on the interface condensa- 
tion rate is enhanced as Peclet number increases or inflow Reynolds number decreases. A correlation based on 
numerical solutions to describe the average condensation Nusselt number, when the effect of buoyancy is not neg- 
ligible, is given by 

Nuc = NUco(l - 0.0062 Ri Pe Re^) (22) 

Equation (22) predict Nu c within ±5 percent of numerical solution for the range of Pr and Re in covered in the 
present study. 

The subcooled liquid inflow introduced upward from the outer portion of the tank turns radially inward to the 
tank centerline. The dimensionless radial velocity at the interface (v*^ for various Richardson number is shown 
in Fig. 5. For Ri = 0, v* s is negative everywhere, implying that all subcooled liquid inflow moves toward the 
tank centerline. As buoyancy effect increases, i.e., Ri > 0 the radial velocity at the interface moving toward the 
tank centerline decreases. If the buoyancy is large enough, surface velocity in the centerline region reverses direc- 
tion and moves outward toward the tank wail as shown in Fig. 5 for Ri = 4.4. This means that there may exist a 
small clockwise vortex region underneath the interface in the centerline region. The corresponding dimensionless 
temperature gradient distribution at the interface is shown in Fig. 6. Except for the outflow region near the center- 
line, for Ri = 0, the higher surface velocity results in higher temperature gradient over the interface area than that 
for Ri = 4,4. Also, for the case of zero Richardson number the flow leaves the tank at a higher temperature in the 
outflow region near the centerline, resulting in a lower interface temperature gradient in the centerline region. On 
the other hand, the clockwise vortex of subcooled liquid underneath the interface for Ri = 4.4 results in more con- 
densation in this small confined region near the tank centerline. However, the average vapor condensation rate at 
the interface is higher for Ri = 0 than that for Ri > 0. Figure 7 shows the comparison between the numerical 
predictions and experimental data of Ref. 7 for Pr - 4.37. The agreement is very good for the inflow Reynolds 
number up to about 150. For Pr = 4.37, the subcooied water inflow with Re in > 150 will yield Richardson num- 
ber smaller than 0.5, resulting in negligible effect of buoyancy on the interface condensation rate. However, 

Fig. 7 shows that the numerical solutions under-predict the condensation rate for all data with Re in > 150. This 
may be due to the fact that for higher inflow Reynolds number the flow near the interface is turbulent. Figure 8 
shows the numerical predictions (Eq. 22) and experimental data (laminar) of Ref. 7 for various Prandtl numbers 
are in good agreement. 


Conclusions 

A numerical solution of the direct contact condensation on a slowly moving horizontal liquid surface is 
obtained. The geometrical configuration and the input conditions used are representative to those of the experi- 
ments of Ref. 7. Numerical predictions of condensation rate for laminar flow are in good agreement with experi- 
mental data. Based on the numerical solutions for the parameter ranges covered in the study, the following 
conclusions are made: 

1. For laminar flow without buoyancy effect (Ri * 0), the average condensation Nusselt number (Nu c ) can be 
predicted by Eq. (21) which shows Nu c is linearly increasing with half power of the Peclet number (Pe 1/f2 ). 

2. The effect of buoyancy on the condensation rate can be described by Eq. (22). The average condensation 
Nusselt number is linearly decreasing with increasing Richardson number. An increase in Peclet number or a 
decrease in flow Reynolds number enhances the effect of buoyancy on the average condensation Nusselt number. 
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Figure 1 .—Physical system and coordinates. 


Figure 2. — Effect of Prandtl number on condensation 
rate for Ri - 0. 
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